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At finite temperature, the hydrodynamic collective modes of superfiuid trapped Fermi gases are 
coupled to the motion of the normal component, which in the BCS limit behaves like a collisionless 
normal Fermi gas. The coupling between the superfiuid and the normal components is treated in the 
framework of a semiclassical transport theory for the quasiparticle-distribution function, combined 
with a hydrodynamic equation for the collective motion of the superfiuid component. We develop a 
numerical test-particle method for solving these equations in the linear response regime. As a first 
application we study the temperature dependence of the collective quadrupole mode of a Fermi gas 
in a spherical trap. The coupling between the superfiuid collective motion and the quasiparticles 
leads to a rather strong damping of the hydrodynamic mode already at very low temperatures. At 
higher temperatures the spectrum has a two-peak structure, the second peak corresponding to the 
quadrupole mode in the normal phase. 



PACS numbers: 03.75.Ss,03.75.Kk,67.40.Bz 



I. INTRODUCTION 

Most of the current experiments involving trapped 
atomic Fermi gases focus on the BEC-BCS crossover. 
By changing the magnetic field around a Feshbach reso- 
nance, the scattering length a of the atoms can be var- 
ied from small positive values through very large val- 
ues near the resonance to small negative values. For 
a > 0, fci?a<Cl (where kp denotes the Fermi momen- 
tum) the system can be considered as a Bose-Einstein 
condensate (BEC) of diatomic molecules. The crossover 
region, kp\a\ > 1, is not yet very well understood from a 
theoretical point of view. Finally, on the other side of the 
resonance, when a < 0, kp\a\ <C 1, the system should be 
in the BCS phase if the temperature is sufficiently low. 
However, the BCS critical temperature T c is extremely 
low, and very soon the magnetic field reaches the point 
where T c becomes smaller than the actual temperature 
T, and the system undergoes the phase transition to the 
normal (non-superfluid) phase. 

One possibility to study the crossover experimentally 
is to measure the properties of certain collective oscilla- 
tions. For example, the radial and axial breathing modes 
of a cigar-shaped trapped Fermi gas have been measured 
over the whole crossover region [3, • In these experi- 
ments one can observe how the frequencies and damping 
rates of the modes change from what one expects for a 
BEC to what one expects for a collisionless normal Fermi 
gas. Assuming that, except in the collisionless normal 
phase, hydrodynamics is valid, the measured frequencies 
can give some information on the equation of state in the 
crossover region. 

However, this schematic picture is not completely ac- 
curate. Since the system is in a trap, there is no sharp 
transition from the superfiuid to the normal phase. This 
can be seen as follows: The BCS critical temperature T c 
depends on the atom density p, and the density depends 
on the position r. In the center of the trap, the density 



p(r) and hence the local critical temperature T c (r) are 
higher than in the outer part of the trap. As a conse- 
quence, for a given temperature, the outer part gets al- 
ready normal at a magnetic field where the inner part is 
still superfiuid. To be more precise, a system in the BCS 
phase at finite temperature behaves effectively like a mix- 
ture of superfiuid and normal components with densities 
p s and p n , respectively, which become p s = p, p n — 
in the limit T — and p s = 0, p n = p in the limit 
T > T c . As a consequence, if < T < T c (r = 0), the 
superfiuid inner part of the trap behaves like a mixture of 
normal and superfiuid components, while only the outer 
part with T c (r) < T is completely normal 

If the collision rate was high enough, also the nor- 
mal component of the gas would behave hydrodynami- 
cally. Such a system could be described by Landau's two- 
fluid hydrodynamics which has been applied to collective 
modes in trapped superfiuid gases at finite temperature 
However, although in the recent experiments the 
transition to the normal phase seemed to occur at a value 
of kp\a\ « 2 [l| (i.e., the BCS phase has not really been 
reached), the system behaved already like a collisionless 
normal Fermi gas. Hence it seems to be clear that the 
normal component cannot be treated in terms of hydro- 
dynamics, but a description in terms of a Vlasov equation 
is required. 

We note that there are other approaches to the de- 
scription of the collective modes at finite temperature. 
In particular, let us mention the quasiparticle random 
phase approximation (QRPA) which can be seen as 

the linearized form of the time-dependent Bogoliubov-de 
Gennes (BdG) equations. However, for practical reasons 
this method is limited to systems with spherical sym- 
metry and numbers of particles up to a few times 10 4 . 
Another disadvantage of this method is that it does not 
allow to include a collision term. 

For the case of clean superconductors, a semiclassi- 
cal transport theory taking the coupling between nor- 
mal and superconducting components into account has 
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been developed by Betbeder-Matibet and Nozieres 
Transport theories of this type have also been used for 
describing the dynamics of superfluid 3 He @, Q. In a 
preceding paper [101 ] , we derived the semiclassical trans- 
port equations for the case of trapped atomic Fermi gases 
and applied them to the quadrupole mode of a gas in a 
spherical trap. We found that the presence of the normal 
component leads to a strong damping of the hydrody- 
namic collective mode. The same mechanism might ex- 
plain the strong damping observed experimentally near 
the transition to the collisionless behavior However, 
in Ref. [Icj we had to replace the gap A(r) by a constant 
in order to find an analytical solution of the transport 
equations. Due to this simplification, which cannot re- 
ally be justified, the damping of the hydrodynamic mode 
at a given temperature was much weaker than that ob- 
tained in QRPA calculations [g|. 



In the present paper we will work out a numeri- 
cal method which allows us to treat the realistic in- 
dependence of the gap. In addition, the method is very 
versatile and allows to treat much more general cases 
than can be solved analytically in the constant-gap ap- 
proximation. The basic idea is to replace the continuous 
phase-space distribution function of the quasiparticles by 
a sum of a finite number of delta functions in phase space, 
called "test particles." In the normal phase, the test- 
particle method is routinely used for solving the Vlasov 
equation, e.g. for simulating heavy-ion collisions in nu- 
clear physics [ll| . It has also been applied to the simula- 
tion of the dynamics of normal trapped atomic Fermi 
gases with collision term [l2| and of Bose-Fermi mix- 
tures [HI]. However, to our knowledge, the test-particle 
method has not yet been used in the context of super- 
fluid systems, and in fact the numerical difficulties are 
quite different from those encountered in the usual ap- 
plications. 



The article is organized as follows. In Sec. [TT1 we give 
a brief summary of the transport equations for the BCS 
phase and their linearization in the case of small devi- 
ations from equilibrium. We also give arguments why 
some terms which appear in the equations can be ne- 
glected. In Sec. EH we introduce the test-particle method 
for the case of small oscillations around equilibrium. We 
describe in detail a number of tricky points we encoun- 
tered during the implementation of the method, in par- 
ticular the calculation of the test-particle trajectories, the 
generation of the test-particle distribution in phase space, 
and the initialization after a delta-like perturbation. In 
Sec. IIVI we present the first results obtained with the 
help of this method, again for the quadrupole mode in a 
spherical system. Finally, in Sec. |Vl we summarize and 
draw our conclusions. 



II. TRANSPORT EQUATIONS FOR THE BCS 
PHASE 

A. Summary of the kinetic equations 

In this subsection we will give a brief summary of 
the kinetic equation approach developed by Betbeder- 
Matibet and Nozieres [3] for the case of clean supercon- 
ductors and adapted to the case of trapped atomic Fermi 
gases in Ref. [10| . In this paper we will only give the fi- 
nal equations. For details of the derivations, see Ref. [Toj 
and the footnote [27l ] . 

We consider a dilute gas of fermionic atoms of mass 
to in two equally populated hyperfine states | and |, 
trapped by an external potential V ex t and interacting 
via an attractive short-range interaction which leads to 
a scattering length a < 0. The corresponding classical 
mean-field hamiltonian (minus the chemical potential fi) 
reads 

h(T,p) = £ + V(T)-n, (1) 

where V denotes the sum of the external and the Hartree 
potential, 

V(r) = V ext (r) + V H artree (r) = V ext (r) + gp(r) . (2) 

In the latter equation, g = 4irh 2 a/m denotes the cou- 
pling constant and p is the density per spin state. The 
Vlasov equation (without collision term) for the phase- 
space distribution function g(r, p) in the normal phase 
can be written in the compact form 

g = {ft, g} , (3) 

where {•, •} denotes the Poisson bracket. One way to 
derive this equation is to perform a Wigner-Kirkwood 
expansion up to order ft of the time-dependent Hartree- 
Fock equation pHIT3|. 

In the superfluid phase the derivation of an analo- 
gous transport equation is much more complicated due 
to the presence of the complex order parameter (gap) 
A(r) whose phase describes the collective motion of the 
Cooper pairs. In addition to the density matrix g, there 
exists now an anomalous density matrix (pairing tensor) 
k. The gap A and the anomalous density are related by 
the gap equation 




The time-dependence of g and re is governed by the time- 
dependent Hartree- Fock-Bogoliubov or BdG equations. 
As in the normal phase, the semiclassical transport the- 
ory can be derived from these equations by performing a 
Wigner-Kirkwood expansion up to order ft. However, it 
turns out that it is necessary to introduce a gauge trans- 
formation with a phase <fi(r) that makes the order param- 
eter A real. This corresponds to a transformation into 
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the local rest frame of the Cooper pairs, which are mov- 
ing with the collective velocity v co;; (r) = — (7i/m)V<^(r). 
The effect of this transformation is to change the gap A, 
the single-particle hamiltonian h, the normal and anoma- 
lous density matrices g and k according to 

A(r) = A(r)e 2 ^ (r) = |A(r)|, (5) 

ft(r 1 p) = fc[r 1 p-W^(r)]-^(r) ! (6) 
e(r,p) = e[r,p-fiV0(r)] J (7) 
£(r,p) = K (r,p)e 2 ^. (8) 

Roughly speaking, the phase <fi determines the dynam- 
ics of the superfluid component of the system, while the 
dynamics of the normal component, consisting of ther- 
mally excited quasiparticles, has to be described sepa- 
rately. The distribution of these quasiparticles, denoted 
by u(r, p), obeys the following equation of motion: 



v = {E,v}. 



(9) 



This equation looks formally very similar to the Vlasov 
equation (J3]), except that the hamiltonian h is replaced 
by the quasiparticle energie E, which is defined as 



E 



h 2 ev + A 2 + hod 



(10) 



Throughout this article, the indices "ew" and "od" de- 
note the time-even and time-odd parts of a phase-space 
function, i.e., the parts which are even and odd in p, 
respectively. The quasiparticle-distribution function v is 
related to the normal and anomalous density matrices in 
the new gauge, g and k, by 



i ; = ~ TTTT-t 1 ~ + V od 



2E e 



Re k = 



2E P 



-{l-2v ev ). 



(11) 
(12) 



The Vlasov-like equation (J9]) has to be complemented 
with an equation of motion for the phase (p. It turns out 
that <fi has to be determined from the continuity equation 

p(r) + V-j(r) = 0, (13) 

where the density p and the current j are given by 



P(r) 



j(r) 



d 3 p 



(2nh) 



(14) 



/^^g(r,p)--p(rW(r). (15) 



B. Linearization around equilibrium 

Let us now assume that the external potential V ex t can 
be written as 



Vext — Vo ex t + V\ ext j 



(16) 



where Vo ex t is time-independent and V± ex t is a small per- 
turbation. The equilibrium quantities (corresponding to 
the potential Vbext) will be marked by an index "0". In 
particular, since in equilibrium the gap can be chosen to 
be real, we have 



= , ho = hn 



An = Ar 



(17) 



The quasiparticle distribution function in equilibrium is 
given by 



M T >P) = f[Eo(r,p)] , 
where f(E) denotes the Fermi function, 

1 



f(E) 



e E/(k B T) + i 



(18) 



(19) 



Our aim is to calculate the small deviations from equi- 
librium induced by the perturbation V\ ex t, which will be 
marked by an index "1" . To that end we linearize the 
transport equation §§§ for the quasiparticle-distribution 
function, 



(20) 



where f'(E ) = df/dE . The deviation of the quasiparti- 
cle energy, E±, which appears on the r.h.s., depends itself 
on V\ through the deviation of the Hartree field, gpi, and 
the deviation of the gap, Ai. Expressing p\ and Ai in 
terms of v\ , we can write Eq. (|20p as 



!>1 - {-Eo,^l} = 

f'(E ) 



( _ v lext + 
(-P-V 1 



Viext + gpiu - h(j)i 



9 A 



A T7 A (Fi eKf +gpi v - h<fii) 



1 + gA 



h A Ai„ h h 2 
Eft gA m Eq 

-h^-(VV )- V^i - ^(VAo) ■ V0i) . (21) 
bo &o ' 

where p\ v and Ai^ are the quasiparticle contributions to 
pi and Ai, 

f d 3 p h (r,p) 

*^~> (23) 

while A(r) is a function which depends only on equilib- 
rium quantities. The explicit expression for the function 
A(r) reads 



mp F {r) 



(24) 



where the local Fermi momentum Pf( v ) is defined as 
usual by p F (r) / (2m) = e F (r) = p, — Vo(r), and the tem- 
perature dependence of A(r) is governed by the function 



(25) 
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with E$ = ^/£, 2 + Aq(t). In the two limiting cases T = 
and T > T c (r), the function ip(r) takes the values and 
1, respectively. As a consequence, A(r) = if T > T c (r). 

In order to determine the phase 0i, we also linearize 
the continuity equation (fT3|) : 



(26) 



p 1 (r) + V-j 1 (r) = 0. 



Again, we express all quantities in terms of equilibrium 
quantities and the unknown quantities v\ and (f>i. Ac- 
cording to Eq. (|15p . the current j ± can be decomposed 
into quasiparticle and superfluid contributions, 



jl(r)=ji„(r)--po(r)V^i(r) 

m 



0. 



where the quasiparticle contribution is given by 



hui T ) 



d 3 p p 



fi(r, P) 



(27) 



(28) 



(note that only the time-odd part of V\ contributes to the 
integral, and V\ & = Q\ d). The time derivative ii\ which 
appears when one writes down the explicit expression for 
pi can be eliminated with the help of Eq. (|2ip. As a 
result, the l.h.s. of the continuity equation becomes 



A 



1 + gA 



2tT 2 H 3 

mpp 



g — V- (poV0i) 



+ <?V-j u 



Ao. 



d 3 p A p 

(2tt/i) 3 £ 2 



(29) 



As noted in Ref. [l(J, the continuity equation is trivially 
satisfied in the normal phase (T > T c ). This becomes 
evident if its l.h.s. is written in the form (|29p . since in 
the normal phase we have Aq = and A = 0. 



C. Identification of important and unimportant 
terms 

Eq. (f2"Tj) is still very complicated. In order to sim- 
plify the problem, let us look more closely at the different 
terms in order to see if some of them are less important 
than others. The basic assumption being that A, ksT c , 
and ksT are much smaller than e F , the distribution func- 
tion is sharply peaked near the Fermi surface. Under this 
condition it is useful to express the distribution function 
v in terms of the variables r, £, and p instead of r and 
p, where 

£ = /Jo(p,r)«VF(r)[|p|-p F (r)], p^, (30) 

and Vp(y) = pp(r)/m. In terms of these variables, v is 
sharply peaked near £ = 0, and the relevant values of £ 
are of the same order of magnitude as A, fcsT c , and ksT. 



If v\ is written as a function of the new variables, the 
Poisson bracket on the l.h.s. of Eq. (|2"Tj) becomes 



{£0,^1} 



A n , ,cVi 
^P-(VAo) — 



Er 



Aq^ 

Eq p F 



Eq p F 

(VAo) 



(Wo) 
dp 



with the short-hand notation 
dv\ dv\ dp 



dp 



d-d p d-d p 



sin' 



1 dv\ dp 

■dp difp d(p p 



dv\ 
dp 

(31) 



(32) 



where "d p and (p p denote the angles characterizing the 
unit vector p = (sin -d p cos tp p , sin "d p sin (p p , cos i? p ) . 

In addition to the assumption A <C cf, our semiclas- 
sical theory requires that all quantities vary slowly in 
space, i.e., on a length scale L which should be larger 
than the coherence length hvp/iit/S.) [l5j]. Then, us- 
ing A - Eq - £ ~ A, V - 1/L, d/d£ ~ 1/A, and 
d/dp ~ 1, all terms in Eq. IpTT]) can be estimated to be 
of the order of magnitude (vp/L)vi, except the last one, 
which is of the order (vp / 'L)(A/ 'ep) v i- Hence, the last 
term of Eq. (f3"Tj) is negligible. 

Let us now distinguish different kinds of contributions 
to v\, depending on whether they are even or odd func- 
tions in £ and p: 



the part of vi which is odd in £ and even in p 
describes, roughly speaking, a change of the Fermi 
momentum, i.e., fluctuations of the density, and 
contributes to p\ u , 



Plu 



mpF 
2ir 2 h 3 



17 / d ^ 



(33) 



fie 



with d£l p — siwdpd'&pdipp, while its contribution to 
Ai u is suppressed by one power of A/e F and can 
be neglected. 

the part of V\ which is even in £ and odd in p 
describes a shift of the Fermi sphere and therefore 
contributes to the current 



Jii/ 



Pi 



(34) 



and also to the other integral in the continuity 
equation ([29]) . 



v\ ee : the part of V\ which is even in £ and in p describes, 
roughly speaking, a local temperature fluctuation 
and leads to a non- vanishing value of Ai (via Ai„), 



Ax 



1 



1 



dSlp f A 
17 / * -E^ 



(35) 



V\ 



while its contribution to p\ v is suppressed by one 
power of A/ef and can be neglected. 

the part of v\ which is odd in £ and odd in p gives 
only a negligible contribution to the current 
(suppressed by one power of A/ tp)- 
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If one neglects the last term in Eq. (|3~Tj) . the Poisson 
bracket in Eq. (f2"Tj) leads only to a coupling between v leo 
and v\ oe and between v\ DO and v\ ee - To be more specific, 
v\ ee and v\ 00 do not contribute to the dynamics of V\ oe 
and V\ eo . Since we are interested in density oscillations 
and currents, which are determined by v\ oe and v\ eo , we 
might wonder if we could disregard completely v\ ee and 
v\ OQ . To that end we have to check that also on the r.h.s. 
of Eq. (|2ip there is no term which couples the undesired 
quantities v lee and v loo to z/ loe or v leo . Actually, on the 
r.h.s. of Eq. (|2ip there is no term containing v\ 00 and 
only one term containing z^i ee , namely the third one, 



space, 



f(E Q ) h a A0A1, 



m 



P2J 

E 



v F f'(E )^p-V(A A 1 ). 



This term clearly contributes to v\ c» but at least to lead- 
ing order in A/ep it does not contribute to i>i eo or v\ oe . 
In the continuity equation (|29[) . v\ ee and v\ QO do not ap- 
pear, i.e., the undesired quantities v\ ee and v\ ao do not 
contribute to the dynamics of <p\ , either. We are there- 
fore allowed to disregard them. 

Now, since we are not interested any more in v\ ee and 
viooi we can remove all the terms on the r.h.s. of Eq. (|21|) 
which contribute only to the dynamics of these uninter- 
esting quantities. As mentioned above, this is the case 
for the third term, Eq. (f3"6")> . which contributes only to 
v\ 00 - The last term on the r.h.s. of Eq. (f2"Tj) , 



m k,Q 



(37) 



can be omitted, too, since it contributes only to v\ e e- 
In conclusion, we are left with a simplified version of 
Eq. (|2Tj) , which reads 



v\ - {E , v\} = 
f'(Eo) 



h ho 
to En 



P 1 + 9A 

A Q (V lext + gpiu - hipx) 
1 + gA 

£jq 



u x-r 

E? P ' V 



(38) 



III. TEST-PARTICLE METHOD 

A. Description of the method 

The aim of the present work is to solve the Vlasov-like 
equation ([9]) for the quasiparticle-distribution function v 
together with the continuity equation (|13p for the phase 
of the order parameter with the help of the test-particle 
method, in analogy to the test-particle method which is 
used to solve the usual Boltzmann equation. The basic 
idea of this method is to replace the continuous distribu- 
tion function v(r, p) by a sum of delta functions in phase 



Kr, p;t) oc S i r ~ RiitMP ~ Pi(*)] , (39) 

i 

corresponding to a finite number of test particles, each 
of which follows the classical equation of motion 



dE(R z ,Pf,t) 



dg(R„p^) 



(40) 



as can be seen by inserting Eq. (|39|) into Eq. ([9]). Note 
that, contrary to the usual test-particle method, our test 
particles here cannot be identified with real particles but 
rather with Bogoliubov quasiparticles. In its general 
form, the test-particle method can be applied to situa- 
tions far from equilibrium. However, here we are only in- 
terested in the linear-response regime, i.e., in the limit of 
small deviations from equilibrium. In this case it is possi- 
ble to formulate the method in such a way that only the 
classical trajectories corresponding to the unperturbed 
system appear. 

To that end, we make the following ansatz for the de- 
viation of the distribution function from equilibrium: 



(41) 



fl(r,p;i) = -y(r,p;t)f'[E (r,p)] . 

Inserting this into the linearized transport equation (|38[) , 
we obtain the following equation of motion for the func- 
tion y: 

y(r, p; t) - {E {r, p), y(r, p; *)} = F(r, p; t) , (42) 
where 



F(r>p;i ) = £ . v *W + gpi„ -^i 

to 1 + gA 

AoP ^ Ag(V lext +gpi„ - h<fii) 
Elm 1 + gA 

+, £(£-*) , *-s£™- v *- (43) 

Denoting by R(r, p; t) and P(r, p; t) the classical trajec- 
tories satisfying the equations of motion 



R 



d£o(R,P) 
dP 



dE (R,P) 
dR 



with the initial conditions 

R(r,p;Q)=r ) P(r,p;0) = p, 
one can easily show that 



(44) 



(45) 



-y[R(r, p; t), P(r, p; t); t] = F[R(r, p; t), P(r, p; t); t] . 

(46) 

Let us now replace the quasiparticle-distribution func- 
tion by N v delta functions in phase space. Since the order 
of magnitude of v\ is dominated by —/'(Eq), it is clear 
that these delta functions should be distributed near the 
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Fermi surface. To be more specific, we choose N v points 
rj, in phase space which are distributed according to 
a probability density which is proportional to —f'(E ), 
in such a way that for arbitrary but sufficiently smooth 
phase-space functions x( r : p) the integral of x(r, p) times 
the function f'[Eo(r, p)] can be approximated by 

r w 3 w 3 Nv 

j ^.Pim^p^-c&^Pi)- (47) 

Note that, if r i5 p^ are distributed in such a way, the 
same is true for Rj(i) = R(r i ,p i ;f), P,-(t) = P(r i ,p i ;i), 
since the quasiparticle energy Ei — £o[Rj(i), P«(i)] is a 
constant of the motion. In particular, defining j/j(t) = 
t/[Rj(i),Pj(i);t] and using Eq. (|47[> . we can approximate 
the integral of an arbitrary function \ times the distri- 
bution function v\ as 



d 3 rd 3 p 
(2nh) 3 



N„ 



X (r, p)^(r, p; t) w C ^ tfiWxP-i (*), P<(*)] • 



i=l 



In other words, we have replaced v\ by 



(48) 



N,. 

Mr,P;t)^cY^y l (t)S[r~R l (t)}6[p-P t (t)}. (49) 
i=l 

According to Eq. (|46| , the equation of motion of the co- 
efficients t/i is reduced to 

y i (t) = F[R i (t),P i (t);t\. (50) 



Above we assumed the function x(r, p) to be suffi- 
ciently smooth. Of course, this causes some trouble if 
we want to calculate local quantities like the density or 
the current. For instance, we obtain 



N„ 



M*) = c£fi(*)^r*[ r -R*(*)] 



E, 



(51) 



where £j(t) = h [Ri(t), P»(i)]. This result makes sense 
only after the delta functions have been averaged over a 
volume containing a sufficiently large number of test par- 
ticles in order to have a reasonable statistics. Supposing 
that this can be done, and supposing that V ex t(r;t) and 
the phase 0i(r;i) are known, we can use the result for 
p\ v in the explicit expression for F in order to obtain 
a system of N v coupled first-order differential equations 
of the form ([ST))) for the coefficients t/,-. This represents 
a tremendous simplification with respect to the original 
partial differential equation (|38[) in seven dimensions (r, 
p, and t). 

However, the phase <f>i(r, t) is not known, but it has to 
be determined from the continuity equation (|26p . This 
is, again, very difficult. Hence, instead of solving the 
continuity equation exactly, we make an ansatz for <f>i and 
determine the parameters by minimizing the violation of 
the continuity equation, 



d 3 r(pi 



V-Jr) 



(52) 



the explicit expression for pi + V • jj being given by 
Eq. (f2"9"| . The idea is to expand <t>\ on an appropriately 
chosen set of orthogonal functions tl> n , 



(53) 



Inserting this ansatz into Eq. (|29[> . we see that the in- 
tegral in Eq. (|52|) depends on x n and x n . At a given 
time, we regard x n and x n as given (e.g., at the mo- 
ment when the perturbation is switched on, we know 
that x n = x n = 0). Hence, in order to have a minimal 
violation of the continuity equation, we have to minimize 
Eq. (|52p with respect to x n by demanding 



d 



J d 3 r( Pl + V-h) 2 = 0. 



(54) 



At this stage it turns out to be convenient to choose the 
basis functions ip n such that they satisfy the orthogonal- 
ity relation 



d 3 r( - — A Vn( r )^m( r ) = Sum ■ (55) 



1 + gA 



Then we obtain the following differential equation for the 
coefficients x n : 



n(t) = a nmX m {t) + d r .. 2 \Vi ex t 

m=l J K ~>~ 9 ) 



— . A f d 3 p A p 
9V-j lK -— V 



A 7 (2tt/i) 3 £ 2 m 
where a is a time-independent matrix, 



^i), (56) 



A 2 ip n (2-K 2 h 3 



m J (1 + gA) 2 V mpp 

(57) 

Using Eq. (|48p and integrating by parts, we can rewrite 
Eq. (|56p in a more convenient form as 



m— 1 z— 1 

where b(t) denotes the matrix 



m V (1 



(l+.gA) 2 £ 2 (l+ 5 A) 2 ;R l(t ) 

(59) 

and the vector v is defined by 



v n = H d r 



1 



j3 A 2 ip n V lext 
~ (1 + gA) 2 



(60) 



Mainly for formal purposes, we note that also the equa- 
tion (|50p for the coefficients yi can be rewritten in matrix 
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notation as 



Vi(t) = y^Acin(t)Xn(t) + d m (t)x n (t)} + fi(t) 



where 



A _ A ^„ 



= ~P<(t) • (V^^-^V: 



E? l+gA/R»(t) 



»g( W')'^ - |w ,). % , 



(62) 

(*)' 
(63) 



, m _ gift) /- JW_ _ Ao v AoVW ' 
M J m ' V 1 + gA Ef 1 + gA Ju z (t) 



(64) 



are not dealing with ordinary particles but with Bogoli- 
ubov quasiparticles, which have some surprising proper- 
ties. For instance, Ei being a constant of the motion and 
Ef — £| + Ag(rj), it is evident that the energy £3 cannot 
be conserved if the gap Ao depends on r. In particular, 
when a test particle with quasiparticle energy Ei reaches 
the surface where Ao(r) = JSj, it is reflected (Andreev re- 
flection) . During this reflection, the momentum stays 
almost constant, but the energy £j changes its sign (i.e., 
a particle is transformed into a hole or vice versa), such 
that the velocity Vj = dEi/dPi = (£j/-Ej)P,-/m is re- 
versed. As a consequence, the quasiparticle is reflected 
into the direction where it came from, which is very sur- 
prising if the incident angle is different from 90°. 

In order to find the test-particle trajectories numeri- 
cally, it does not seem very efficient to start directly from 
Eqs. (|44p , since a small numerical error in the momentum 
of the order of 8P/P ~ A/ep would immediately lead to 
a completely wrong behavior. It is therefore advanta- 
geous to make use of the variable whose equation of 
motion reads 



and 



9ij(t) = -gC- 



&(*)Pj(t) / V jjr-Rjft)j 
Ej m V 1 + gA 

KM 



A — Apflr- 
E? h 



Ri(t) 



(65) 



In the latter equation, S denotes a kind of "smeared" 
delta function which accounts for the averaging men- 
tioned below Eq. (f5Tj) . 

However, as mentioned above, Eqs. (jUT]) and will 
be used for formal purposes only. In practice, it is much 
faster to calculate pi v (r) after each time step on a dis- 
crete mesh, and to interpolate the stored values when 
performing the next time step for the coefficients y^. For 
the calculation of p\ v (r) , we replace the delta function in 
Eq. (|51|) by a Gaussian having a width d p . 

In summary, the coupled system of partial differential 
equations, namely the transport equation for the distri- 
bution function v\ and the continuity equation for the 
phase 4>\ [Eqs. (j3"5)) and (|26p ], has been replaced by a 
coupled system of ordinary linear differential equations 
for the coefficients yi and x n [Eqs. (|6T|) and (I58)) ]. which 
can formally be written as 



d 



f 1 

a b(t) 
K d(t) c(t) g(t) 




(66) 



B. Trajectories of the test particles 

In practice, the solution of the classical equations of 
motion for the test particles, Eqs. (144|) . faces us with 
some unusual features which are not present with the 
usual Newtonian equations of motion. Note that we 
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A (Rj) P a 

Ei m 



VA (Pm) . 



(67) 



Solving this equation together with the equations for R; 
and P^ we can correct Pi after each time step according 
to 



(68) 



In practice, the variable also allows us to introduce 
a very reliable method for determining the step size. Let 
us denote by £ z ' the result we obtain after one time step 
of size St, and by the result we obtain after two time 
steps of size St/ 2 each. Then the quantity St\£' — £"| is 
a measure for the numerical error and can be used for 
adapting the step size St to the situation. It turns out 
that the step size has to become very small only during 
Andreev reflection. 

Now let us give some examples for typical test-particle 
trajectories. For that purpose, let us restrict ourselves to 
the most simple case which is a spherical harmonic trap, 



V 0ext (r) = ^mtfr 2 



(69) 



This potential defines the so-called trap units, i.e., ener- 
gies are measured in units of Ml, temperatures in units 
of Ml/kg, lengths in units of Iho = y/h/ (mil), etc. In 
this example, due to spherical symmetry, not only the 
quasiparticle energy E, but also the angular momentum 
L = r x p of a test particle is a constant of the motion. 

Within the local-density approximation (LDA) (ltl 
f]~7| . the density po(r) has its maximum at the center 
of the trap and vanishes approximately (except for very 
small temperature effects) at the Thomas-Fermi radius 
Rtf — y^2p/(mfl 2 ). The gap Ag(r) has its maximum 
at the center of the trap, too, and goes to zero at some 
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critical radius R c which is temperature dependent and 
determined by the equation T = T C (R C ). In order to 
avoid numerical problems arising from the infinite deriva- 
tive of Ao(r) at r = R c , we convolute the LDA result for 
Ao(r) with a Gaussian of width c?a- In fact, this is more 
realistic than the LDA result since the exact solution of 
the BdG equations also leads to a gap Ao(r) which has 
an exponential tail [l7|, |l8|, [l9( . As parameters we choose 
[i = 32 Ml, g = -h 2 l ho /m, and T = 1.4 Ml/k B . The cor- 
responding number of atoms in the trap is approximately 
1.7 x 10 4 . For these parameters quantum mechanical 
(BdG, QRPA) results are available for comparison. The 
width c2a is chosen such as to optimize the agreement 
with the BdG gap, which for the present parameters is 
achieved with c?a = I ho- 

In Fig. [T] we show the corresponding gap Aq(t) as a 
function of the distance r from the center of the trap. 
From this figure it is evident that due to the condi- 
tion E > Ao(r), the relevant quasiparticles (having E < 
ksT = 1.4 Ml) are excluded from the region r < 4/; lo . In 
addition to the gap, we display the potential Vo(r) — fi, 
since the motion of a quasiparticle with given energy E 
and angular momentum L is also limited by the condi- 
tion ^JE 2 - Ag(r) - L 2 /2mr 2 > V (r) - \x. It has been 
shown that also within the fully quantum-mechanical 
BdG theory the lowest-lying quasiparticle states are lo- 
calized in this region [20(. In our example, the motion 
of the relevant quasiparticles is restricted to the region 
4<r/J ho <8. 

Most of these quasiparticles will undergo Andreev re- 
flection. Their trajectories are approximately described 
by an ellipse which is cut at the points where Ao(r) = E. 
If E <C £p, the quasiparticle will move hence and forth 
on the same partial ellipse. Such trajectories with E = 
0.1 Ml and 0.4 M7 are shown in the left panel of Fig. [5] 
However, if the quasiparticle energy is higher, the change 
in energy from ^ w £ to f w - E (or vice versa) during 
the Andreev reflection results in a change of momentum 
which is no more negligible. Then, due to angular mo- 
mentum conservation, the angle of reflection is slightly 
different from the angle of incidence, and the whole tra- 
jectory is precessing. An example for such a trajectory 
with E = 0.7 Ml is also shown in the left panel of Fig. [2] 
A completely different picture arises if the initial condi- 
tions are such that the quasiparticle does never reach the 
point where A(r) = E. Then the trajectory is just a pre- 
cessing, slightly deformed ellipse, as shown in the right 
panel of Fig. [2] for the case of a trajectory with E = Ml. 
There is a striking analogy between these trajectories and 
the "glancing" orbits discussed, e.g., in Ref. [2l| in the 
context of a superconducting cylinder which is coated by 
a normal- metal layer. 



C. Distribution of test particles in phase space 

In Sec. IIII Al we supposed that one can generate a dis- 
tribution of points Yi, p i in phase space such that Eq. (14"?| 
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r Oho) 

FIG. 1: Gap Ao(r) (solid line) and potential Vo(r)—fi (dashed 
line) for the case of a spherical trap with frequency f2, chem- 
ical potential fj, — 32 Ml, coupling constant g — —h 2 lho/m, 
and temperature ksT =1.4 HQ. Ao and Vq — fi are in units of 
HQ, r is in units of the oscillator length i; lo . Roughly speak- 
ing, these two curves determine the classically allowed region 
for a quasiparticle with given energy E. 
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FIG. 2: Four examples of quasiparticle trajectories in a trap 
with parameters given below Fig. [TJ The three trajecto- 
ries shown in the left panel belong to quasiparticles with 
E — 0.1 Ml, 0.4 M7, and 0.7 Ml, respectively. The trajectory 
displayed in the right panel belongs to a quasiparticle with 
E = Ml. 

is approximately satisfied for sufficiently smooth func- 
tions x( r : p)- I n practice, this distribution is obtained in 
two steps. First we generate the coordinates r^, and in a 
second step the momenta p.^. 

The mean density of test particles at a certain point r 
is given by 

n(r)=5^(ri-r), (70) 

i=l 

where S denotes a smeared delta function in order to 
account for the averaging. Using Eq. (|4"T|) . we conclude 

The algorithm for the generation of the coordinates is 
now very simple. First we look for the maximum w rnax 
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of the function w(r). Defining P(r) — w(r)/w max , we 
obtain a function whose values lie between and f . Then 
we generate uniformly distributed random points in 
a volume which contains the whole system, and retain 
each point with the probability P(?k), until the desired 
number of points, N u , is reached. 

The formula (f7T|) for the test-particle density n(r) can 
also be used for the determination of the normalization 
constant C. Integrating n(r) over space, we must recover 
the total number of test particles. This implies 



C 



lkJ d3rW{r) - 



(72) 



Now we turn to the distribution of the momenta p L . It 
is evident that the angular distribution of the momenta 
is isotropic, i.e., the interesting part of the problem is the 
distribution of the absolute values, pt — jpj, which is, of 
course, directly related to the distribution of the energies 
£j. Let us define the mean number of test particles per 
energy and volume 



n(r,£)=^S( ri -r)£(&-£). 

i=l 

Again, with the help of Eq. (j4"Tj) , this becomes 



(73) 



(74) 



with p£ = ^/2m£ + pp(r), i.e., for given spatial coordi- 
nates r, the probability density for finding a particle at 
energy £ is proportional to —p^f'(E^). Such a distribu- 
tion can be generated in the following way. Starting from 
random numbers Zk which are uniformly distributed in 
the interval (0, 1), it is straight-forward to show that the 
energies 



&=Tln 



Zk 



(75) 



are distributed according to the probability density 
— /'(£)■ It is evident that negative energies with £ < 
— cf(i") have to be removed. Furthermore, it is prefer- 
able to cut the distribution at energies which lie too far 
away from the Fermi surface, e.g., |£| > 15T (the proba- 
bility that this happens is less than 10 -6 ). The momenta 
P£ are thus limited by p ma x = y/S0mT + p F {r), and the 
function defined by P(£) = p$f'(E e )/p max f'(£) cannot 
become greater than 1 and can serve as a probability. If 
we retain each energy generated according to Eq. (|75j) 
with the probability P(£fc), the remaining energies are 
distributed according to the desired distribution. 

In order to give an illustration for the resulting dis- 
tribution of test particles, we show in Fig. [3] the radial 
distribution of N v = 10 5 test particles in a trap with the 
same parameters as in Fig. [TJ In agreement with what we 
discussed in the preceding subsection, we see that the test 
particles are mainly located in the region 4 < r/lf, < 8, 
corresponding to the region where the system is mainly 



r dho) 

FIG. 3: Radial distribution of 10 test particles in a trap 
with parameters given below Fig. [1] counted in bins of size 
Sr = 0.1 Iho- For comparison, the dotted curve represents the 
ideal distribution according to Eq. (J7TJ). 



normal fluid. Due to the angular average the statisti- 
cal fluctuations around the ideal distribution, Eq. (I7ip . 
which is represented by the dotted line, are very small. 
We verified that, apart from the statistical fluctuations, 
our test-particle distribution stays constant, which is a 
good numerical test of both the initial test-particle dis- 
tribution and of the test-particle trajectories. 



D. Initial condition 

In the linear response regime, as the name implies, the 
response to a time-dependent perturbation of the form 
Vi(r;t) = Vi(r)/(£), with an arbitrary time dependence 
fit), can be obtained as convolution of f(t) with the 
response to a delta function in time. It is therefore suffi- 
cient to study perturbations of the form 



V lexf (r;t) = V 1 (r)S(t). 



(76) 



We thus set the inhomogeneous terms in Eq. (f6"6"]) to 
v(t) = v5(t) and f(t) — fS(t), respectively, v and / being 
defined analogously to Eqs. (|6"0|) and (|64p but with Vi e xt 
replaced by V\ ■ 

Assuming that the system was in equilibrium before 
this perturbation, we may ask the question: What are 
the values of the coefficients j/j and x n immediately af- 
ter the perturbation, i.e., at infinitesimally small t > 0? 
This question can be answered exactly, since during the 
infinitesimal period where the perturbation is active, the 
matrix in Eq. (|66[) can be regarded as time-independent. 
Integrating Eq. (|66p over time from — to to to, we obtain 
in the limit to — > 



(77) 




y(*o)> 




Let us now assume that the function V\ lies in the 
space spanned by the functions ip n . Then it is evident 
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that the corresponding linear combination is given by the 
coefficients v„, i.e., 



JV,„ 



(78) 



Note that the functions ip n do not necessarily have to 
have this property. For example, we could define a ba- 
sis of functions satisfying the orthogonality relation ((55) 
and maybe even a suitably defined completeness relation 
if N<p — > oo, but which all vanish identically outside the 
superfluid region, i.e., in the region where A = (and 
A = 0). Eq. I|78p would then be satisfied inside the super- 
fluid region, but not outside. Hence, it is an additional 
requirement for the choice of the functions ip n . Combin- 
ing Eqs. E7J and EE), we find 



lim0i(r;to) = ivi(r) 

CO - fi 



(79) 



Eq. (|78| also leads to a simplification of the initial value 
of the coefficients yi and the quasiparticle-distribution 
function. Using the explicit expressions for the matrix c 
and the vector / [Eqs. ((62) and (164) with Vy e xt replaced 
by Vl], we obtain from the third line of Eq. 1(77) 



lim J/i(*o) = } Ci n V n + fi 
tn—>0 z — * 



Ei.f. 

m V 



h J2n=l ^n^n ~ Vl 

1 + gA 



~Ef V TT^I Jr.' (80) 

As a consequence, if Eq. ((75) is satisfied, the initial values 
of the coefficients yi vanish, which implies 



lim z/i(r,p;f ) = 0. 

r — *0 



(81) 



In fact, the simple result of this subsection, which is 
summarized in Eqs. (|79[) and Ij81). could have been an- 
ticipated without any calculation. The effect of a pertur- 
bation of the form ((75) is to give a particle at position r 
a kick 



5p = - dtVV lext (v;t) = -Wi(r) 



(82) 



Since this kick does not depend on the momentum of the 
particle, the local Fermi sphere is shifted as a whole, there 
is no change in density and no Fermi surface deformation. 
Within the present theoretical framework, Cooper pairs 
are not broken either, they just acquire a center of mass 
momentum. Thus, the distribution function in the local 
rest frame stays unchanged [y\ = 0), and the collective 
velocity is given by v c0 « = — (/l/m)V^i = — (l/m)VVi. 

Note, however, that in reality a perturbation which 
has the form of a short pulse would lead to much more 



complicated effects (e.g., pair breaking). Since our semi- 
classical description requires that the time dependence 
of the perturbation is slow, our formal result for a delta- 
like excitation becomes physically meaningful only after 
it has been convoluted with a function f(t) which varies 
slowly in time. In other words, we can only calculate the 
low-frequency part of the response function. 



IV. FIRST RESULTS 

In this section we will discuss first numerical results 
which have been obtained using the test-particle method. 
Our intention here is to see whether this method is in 
principle capable to describe the most important fea- 
tures of collective excitations in superfluid trapped Fermi 
gases. To that end, we will study the quadrupole excita- 
tion of a spherical system, which is excited by 



Vi(r) = amft(2r 2 z 



;) 



(83) 



(the factor mf2 has been introduced in order to make the 
coefficient a dimensionless) . 

For practical purposes, we will make an additional ap- 
proximation: We will restrict our ansatz for the phase, 
Eq. (|53p . to only one or two functions ip n . It is clear from 
rotational symmetry that in the case of a quadrupole ex- 
citation of the form (|83[) the most general form the phase 
can have is 

1 (r) = $(r)[2rJ-^-^] ) (84) 

such that that the functions ipn can be written as 

^ n {v)=^ n {r)[2rl-rl-rl\. (85) 

It is known from superfluid hydrodynamics that at zero 
temperature the velocity field is essentially linear in the 
coordinates, i.e., the function $(r) is almost constant. As 
a first guess we will assume that this is still true at non- 
zero temperature, and hence we will take only one single 
function = 1) in the ansatz ((53") for the phase, \&i = 
const. The proportionality constant will be determined 
from the normalization condition ((55). 

Such a restricted ansatz means of course that the con- 
tinuity equation will not be exactly satisfied in the super- 
fluid region (remember that outside the superfluid region 
the phase has no effect whatsoever). We will therefore 
improve this initial ansatz by including a second function 
(N0 = 2) which allows to modulate <E>(r) in the superfluid 
region. 

The first idea one might have is to use for ^> n (r) poly- 
nomials in r 2 and to orthogonalizc the resulting functions 
ip n . However, it turns out that this leads to numerical 
instabilities due to the fast growing of the resulting poly- 
nomials outside the superfluid region. Let us explain this 
effect in some more detail. As seen from the transport 
equation for the quasiparticle-distribution function, the 
phase 4>\ outside the superfluid region enters directly the 
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r dho) 

FIG. 4: Density profile po(r) (dashed line) and gap Ao(r) 
(solid lines) in a spherical harmonic trap containing 1.7 x 10 4 
atoms (fj, = 32 Ml, interaction strength g = —h 2 l) lo /m). The 
gap is displayed for three different temperatures, T/T c = 0.2, 
0.4, and 0.6, while the density profile is practically indepen- 
dent of T. 



dynamics of v\ . Although the net effect of the phase and 
of the quasiparticle motion should be independent of the 
choice of <j>\ outside the superfluid region, each of these 
contributions depends on this choice. If <f>\ changes too 
rapidly, the numerical solution of the equation of mo- 
tion for the coefficients yi becomes less accurate and the 
cancellation of the two effects does not work any more. 

We therefore have to look for functions ^> n which are 
linearly independent inside the superfluid region, but 
which do not grow outside. Here we will choose the func- 
tions ^i(r) = 1 and ^(r*) = [1 — l -p{ r )] 2 , where <p(r) is 
the function defined in Eq. - The function ^(r) has 
its maximum in the center of the trap and goes smoothly 
to zero at the boundary of the superfluid region. From 
\Pi and ^2 the functions *i and iff 2 are determined ac- 
cording to the orthogonality condition (|55[) with the help 
of the Gram-Schmidt orthogonalization method. As we 
will see, the results obtained with = 1 and = 2 
are very similar and we therefore claim that they would 
not change a lot if we included additional functions. 

Let us now present the results. As in the examples 
shown in the preceding section, we consider a spherical 
harmonic trap with p = 32M7, containing approximately 
1.7 x 10 4 atoms. The resulting density profile po(r) is 
shown in Fig. Q] as the dashed line. The critical temper- 
ature within LDA is T c = T c (r = 0) » 3.9Ml/k B - As 
before, the LDA result for the gap Arj(r) is convoluted 
with a Gaussian having a width c?a = Iho- We will study 
the quadrupole mode for three different temperatures, 
T/T c = 0.2, 0.4, and 0.6. The equilibrium gap A (r) for 
these three temperatures is also displayed in Fig. [4] 

After the system is excited, its shape will oscillate. A 
measure for this quadrupole deformation is the ratio 



(27 



(r 2 )o 



(86) 



where (r 2 )o denotes the mean square radius in equilib- 



rium, which in the present case has the value (r 2 )o ~ 
23Z^ Q . In the linear response, the quadrupole deforma- 
tion is of course proportional to the strength of the per- 
turbation, and we therefore divide our results by this 
strength [denoted a in Eq. (|83[) ]. In our simulation we 
use N„ = 10 5 test particles, the width of the Gaussians 
used for smearing p\ v (see Sec. lIII A"| is set to d p = Iho- In 
Fig. O we display the time dependence of the quadrupole 
deformation after the perturbation for the three tempera- 
tures mentioned before. The corresponding spectra, ob- 
tained by Fourier transformation, are shown in Fig. [5] 
The results for the two cases — 1 and N$ — 2 are 
displayed as dashed and solid curves, respectively. In all 
cases the two curves are in reasonable agreement, such 
that we can say that the use of N$ = 2 independent 
functions in the ansatz for the phase is sufficient. 

We see that the temperature dependence of the spec- 
trum is highly non-trivial. At low temperatures (up- 
per panel of Fig. [S]), we see essentially the hydrody- 
namic quadrupole mode, which at zero temperature lies 
at u> = V2n HI H3, H|| and which is now damped 
as a consequence of its coupling to the normal compo- 
nent. At higher temperatures (middle of Fig. [6]), a sec- 
ond peak builds up in the spectrum, corresponding to 
the quadrupole mode in the normal phase, which lies 
slightly above uj = 2fi (2f| (for thepresent set of param- 
eters, its frequency is uj w 2.2 O 6]). As the tempera- 
ture approaches T c (lower panel of Fig. [6]), the strength 
contained in this second peak increases, while the hydro- 
dynamic mode, whose frequency is slightly shifted down- 
wards, disappears. These findings are in good agreement 
with quantum mechanical QRPA calculations 6]. 

We note that the damping width of the hydrodynamic 
mode at low temperature is now comparable with that 
found within the QRPA and much stronger than that 
found in our previous work [lfj, where we replaced the 
gap Ao(r) by a constant. The reason is in fact very sim- 
ple: With a constant gap, the fraction p n / po of the nor- 
mal component is independent of r, whereas in the case 
of an r-dependent gap the normal component in the outer 
part of the system is already important at very low tem- 
peratures [|[. 

As the temperature T approaches T c , the quadrupole 
mode of the normal phase (that at uj = 2.2 Q) becomes 
undamped, as it is the case within the QRPA. However, 
even though collisions are strongly suppressed at these 
low temperatures, it should be kept in mind that the 
collision term, which is neglected in the present work, is 
non-zero and its inclusion would lead to a finite lifetime 
of this oscillation, too. 

Finally, let us compare our semiclassical results more 
quantitatively to quantum mechanical QRPA results. In 
Fig. [7] we show the QRPA result of Ref. H for the 
quadrupole excitation spectrum (dotted line) together 
with the semiclassical result we obtain with the same 
parameters (solid line). As one can see, the total nor- 
malization and the relative weights of the two peaks are 
in reasonable agreement. Also the widths of the QRPA 
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FIG. 5: Time dependence of the quadrupole deformation after 
a delta-like perturbation at t = 0. The parameters are the 
same as in Fig. [4] the three panels correspond, from top to 
bottom, to T/T c = 0.2, 0.4, and 0.6. The dashed and solid 
lines correspond to N$ = 1 and N$ = 2, respectively. 



peaks are well reproduced by the semiclassical calcula- 
tion. The main differences are that within the semiclas- 
sical calculation the two peaks lie a bit too high and that 
they are not as well separated as within the QRPA. A 
comparison of the two curves for iV^ = 1 and N$ = 2 in 
the T/T c = 0.4 case shown in the middle of Fig. [6J whose 
parameters are quite close to those of Fig. [JJ suggests 
that the latter effect might be partly due to the restricted 
ansatz for the phase. However, as one can deduce from 
the irregular structure of the QRPA spectrum, even in a 
system with 32000 atoms shell effects, i.e., effects which 
depend on the discrete single-particle spectrum, are still 
quite pronounced. It is clear that such effects cannot 
be reproduced within a semiclassical calculation. In this 
sense the agreement between the two spectra is very sat- 
isfactory, in particular since one can assume that the shell 
effects decrease with increasing number of particles. 
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V. CONCLUSIONS 

In this paper, we developed a numerical test-particle 
method for solving the semiclassical transport equations 
for an ultracold trapped Fermi gas in the BCS phase in 
the collisionless limit. These transport equations take 
into account the coupling between the dynamics of the 
Cooper pairs (superfluid component) and the thermally 
excited Bogoliubov quasiparticles (normal component). 
We developed the method for the case of small deviations 
from equilibrium, so that the test-particle trajectories 
can be calculated in the equilibrium state. Since the test- 
particles describe Bogoliubov quasiparticles rather than 
real particles, the trajectories have very unusual proper- 
ties compared with the trajectories one has to deal with 
when applying the test-particle method to the normal 
Vlasov equation. Our test particles can have the char- 
acter of particles as well as holes, depending on whether 
their energy £ is positive or negative, and they can also be 
transformed from the one into the other if they hit the re- 
gion where the gap A becomes larger than their quasipar- 
ticle energy E (Andreev reflection). Another complica- 
tion as compared with the normal Vlasov equation is that 
the dynamics of the quasiparticles is coupled to the col- 
lective motion of the superfluid component, which is de- 
scribed by the phase <f> of the order parameter. This phase 
has to be determined simultaneously with the evolution 
of the quasiparticle-distribution function by solving the 
continuity equation. In the present work, we make an 
ansatz for <p with time-dependent coefficients, leading to 
an approximate solution of the continuity equation. 

As a first application, we calculated the response of 
a gas trapped in a spherical trap to a delta-like pertur- 
bation of quadrupole form. After this perturbation, the 
shape of the gas shows a damped oscillation. At low 
temperatures, this oscillation is just the hydrodynamic 
quadrupole mode which is damped by its coupling to 
the normal component. With increasing temperatures, 
the extension of the normal component increases, and, 
as a consequence, the normal component can perform its 
own quadrupole oscillation. Since the frequency of the 
quadrupole mode in the normal collisionless Fermi gas is 



higher than that of the hydrodynamic mode, this leads 
to a two-peak structure in the response function. As 
the temperature approaches T c , the strength of the hy- 
drodynamic mode disappears and only the normal mode 
survives. 

The next step will be to apply the method presented 
here to more realistic cases, namely to the axial and ra- 
dial breathing modes of a gas in a cigar-shaped trap con- 
taining a larger number of particles. In fact, the deforma- 
tion and the large particle number do not pose a big prob- 
lem, which is one of the main advantages of the present 
method as compared with quantum mechanical QRPA 
calculations. Another possible application of the method 
is to study the dynamics of a vortex, where already the 
equilibrium situation is characterized by a non- vanishing 
phase of the order parameter. 

However, there are still a number of unsolved problems 
and possible improvements of the method. First of all, 
the collision term [9] should be included, which is an ad- 
ditional source of damping of the collective oscillations. 
As mentioned in the introduction, the possibility to in- 
clude collisions is an important advantage of the present 
method as compared with the QRPA, where collision ef- 
fects cannot be taken into account since this would neces- 
sitate to include four-quasiparticle excitations. Second, 
from a fundamental point of view, the fact that the conti- 
nuity equation is only approximately fulfilled is of course 
unsatisfactory and one should think about another nu- 
merical method for solving the continuity equation. Fi- 
nally, one might ask the question how the present theory 
can be extended to the strongly interacting regime. Un- 
fortunately, this question is up to now completely open, 
since in this regime thermal fluctuations of the order pa- 
rameter, which are not contained in the BdG equations, 
play a crucial role (see, e.g., Ref. [26j). 

Acknowledgments 

The author wishes to thank P. Schuck for numer- 
ous fruitful discussions and the critical reading of the 
manuscript. 



[1] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, 
C. Chin, J. Hecker Denschlag, and R. Grimm, Phys. 
Rev. Lett. 92, 203201 (2004); A. Altmeyer, S. Riedl, 
C. Kohstall, M. Wright, R. Geursen, M. Barten- 
stein, C. Chin, J. Hecker Denschlag, and R. Grimm, 
|cond-mat/0609390| (2006); A. Altmeyer, S. Riedl, 
C. Kohstall, M.J. Wright, J. Hecker Denschlag, and 
R. Grimm, e-print |cond-mat /061 1285] (2006) . 

[2] J. Kinast, S.L. Hemmer, M.E. Gehm, A. Turlapov, and 
J.E. Thomas, Phys. Rev. Lett. 92, 150402 (2004); J. Ki- 
nast, A. Turlapov, and J.E. Thomas, Phys. Rev. A 70, 
05140RR) (2004). 

[3] M. Urban, Phys. Rev. A 71, 033611 (2005). 



[4] E. Taylor and A. Griffin, Phys. Rev. A 72, 053630 (2005). 
[5] GM. Bruun and BR. Mottelson, Phys. Rev. Lett. 87, 

270403 (2001); GM. Bruun, Phys. Rev. Lett. 89, 263002 

(2002). 

[6] M. Grasso, E. Khan, and M. Urban, Phys. Rev. A 72, 
043617 (2005). 

[7] O. Betbeder-Matibet and P. Nozieres, Ann. Phys. (N.Y.) 

51, 392 (1969). 
[8] J.W. Serene and D. Rainer, Phys. Rep. 101, 221 (1983). 
[9] D. Vollhardt and P. Woffle, The Superfluid Phases of 
Helium 3 (Taylor & Francis, London, 1990) 
[10] M. Urban and P. Schuck, Phys. Rev. A 73, 013621 (2006). 
[11] GF. Bertsch and S. Das Gupta, Phys. Rep. 160, 190 



14 



(1988). 

[12] F. Toschi, P. Vignolo, S. Sued, and M.P. Tosi, Phys. Rev. 
A 67, 041605(R) (2003); F. Toschi, P. Capuzzi, S. Succi, 
P. Vignolo, and M.P. Tosi, J. Phys. B 37, S91 (2004). 

[13] T. Maruyama, H. Yabu, T. Suzuki, Phys. Rev. A 72, 
013609 (2005). 

[14] P. Ring and P. Schuck, The Nuclear Many-Body Problem 
(Springer- Ver lag, Berlin 1980), 

[15] A.L. Fetter and J.D. Walecka, Quantum Theory of Many- 
Particle Systems, (McGraw-Hill, New York 1971). 

[16] M. Houbiers, R. Ferwerda, H.T.C. Stoof, W.I. McAlexan- 
der, C.A. Sackett, and R.G. Hulet, Phys. Rev. A 56, 4864 

(1997) . 

[17] M. Grasso and M. Urban, Phys. Rev. A 68, 033610 
(2003). 

[18] MA. Baranov and D.S. Petrov, Phys. Rev. A 58, R801 

(1998) . 

[19] G. Bruun, Y. Castin, R. Dum, and K. Burnett, Eur. 

Phys. J. D 7, 433 (1999). 
[20] G.M. Bruun and H. Heiselberg, Phys. Rev. A 65, 053407 



(2002) . 

[21] C. Bruder and Y. Imry, Phys. Rev. Lett. 80, 5782 (1998). 
[22] G.M. Bruun and C.W. Clark, Phys. Rev. Lett. 83, 5415 
(1999). 

[23] M.A. Baranov and D.S. Petrov, Phys. Rev. A 62, 

041601 (R) (2000). 
[24] M. Cozzini and S. Stringari, Phys. Rev. Lett. 91, 070401 

(2003) . 

[25] L. Vichi and S. Stringari, Phys. Rev. A 60, 4734 (1999); 
C. Menotti, P. Pedri, and S. Stringari, Phys. Rev. Lett. 
89, 250402 (2002). 

[26] A. Perali, P. Pieri, L. Pisani, and G.C. Strinati, Phys. 
Rev. Lett. 92, 220404 (2004). 

[27] There are some typos in Ref. [Tpj ] which we wish to cor- 
rect here: In Eq. (33), A should be replaced by A. In 
Eq. (34b), q ev should be replaced by v m . In Eq. (61), 
Ai„(r) should be replaced by — Ai„(r). In addition, we 
point out that Ai = if Ao = 0, which is not evident 
from Eq. (61) but can be derived from Eq. (54). 



